One-dimensional plasmons confined in bilayer graphene p-n junctions 
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Gapless spectrum of graphene allows easy spatial separation of electrons and holes with an external 
in-plane electric field. Guided collective plasmon modes can propagate along the separation line, 
with the amplitude decaying with the distance to it. Their spectrum and direction of propagation 
can be controlled with the strength and direction of in-plane electric field. In the case of a bilayer 
graphene additional control is possible by the perpendicular electric field that opens a gap in the 
band spectrum of electrons. We investigate guided plasmon spectra in bilayer p-n junctions using 
hydrodynamics of charged electron liquid. 

PACS numbers: 73.23.-b, 72.30.+q 
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I. INTRODUCTION 

Atomically thick graphene^ allows control of its electri- 
cal properties in a variety of ways and thus is of great po- 
tential for nanoelectronics^ and optoelectronics^^ One 
of the major advantages lies with the gapless nature of 
its electron band spectrum, which means that relatively 
moderate fields are required to induce desirable changes 
in the density of electrons and holes. This is in a sharp 
contrast to conventional nanoplasmonics of metal par- 
ticles whose properties are fixed for once at synthesis. 
Plasmons in graphene, due to much higher tunability 
of the latter, acquire new properties. While the plas- 
mon spectrum in a homogeneous graphene sheet sub- 
ject to finite temperature^ or doping^ follows the con- 
ventional two-dimensional (2D) form^ w 2 (q) cx q , spa- 
tial separation of electrons and holes leads to novel ex- 
citations. In particular, one-dimensional plasmons can 
propagate along a graphene p-n junction^ with the spec- 
trum uj 2 (q) cx {E^q) 1 / 2 , which depends on the strength 
of electric field Eq that separates electrons and holes. 
Thus both the direction of propagation and velocity (for 
a given wavelength) of these guided plasmons can po- 
tentially be controlled by changing the orientation and 
magnitude of the electric field that creates the p-n junc- 
tion. 

Bilayer graphene^ is potentially of even greater signifi- 
cance for applications because of the possibility of open- 
ing a bandgap via external gating^ and thus transi- 
tioning from the metallic to the insulating state. Cor- 
respondingly, plasmon excitations in bilayer graphene 
are a subject of active research*^— However, to date, 
only 2D plasmons in homogeneous bilayers have been 
studied. In the present paper we consider guided one- 
dimensional plasmons propagating along p-n junctions 
in bilayer graphene, in an extension of the work in 
Ref. H|. Bilayer p-n junctions have recently been reported 
experimentally^ and their transport properties are at- 
tracting increased theoretical interest J£r— Here we ad- 
dress their collective modes. 

Our objective is to study the behavior of one- 
dimensional plasmons propagating in a bilayer graphene 
and the dependence of the plasmon spectrum on the elec- 



tric field Eq. In this work we consider two cases; the 
first case is plasmon propagation in a bilayer graphene 
flake under the effect of only the in-plane electric field 
E (gapless case), discussed in Sees. II and III. In the 
second case, discussed in Sec. IV, there is an additional 
electric field perpendicular to the plane of graphene that 
opens a gap in the band electron spectrum. In both 
cases the system is a flake of width 2d (in the x direc- 
tion) and infinite length along the y axis, as shown in 
Fig. 1(a). The in-plane electric field E is applied along 
the x axis. Due to the gapless nature of the spectrum 
of bilayer graphene, electric field induces charge separa- 
tion and produces charge density po(x) across the width 
of the flake such that one half of the flake is positively 
charged and the other half is negatively charged. Fluctu- 
ations of charge density Sp propagate on top of the equi- 
librium density in the form of plasmon waves. Because 
of the nonuniform profile of po(x) plasmons are localized 
near the neutrality line x = and behave as quasi-one- 
dimensional excitations. We utilize the Thomas-Fermi 
equation to find the equilibrium charge density pa and 
then use hydrodynamic equations to describe the dynam- 
ics of charge oscillations Sp(r,t). It reduces to solving 
an integrodifferential eigenvalue problem for the plasmon 
frequencies. For the ungated case in the limit of short 
wavelengths -C d even eigenfrequencies have higher val- 
ues than the odd ones. Interestingly, at large wavelengths 
d this order is changed and the lowest cigenmode is 
an even one. In all cases the plasmon frequencies are 
proportional to y/E~o~. In the gated case, the perpendicu- 
lar electric field and the ensuing energy gap leads to the 
appearance of a neutral strip at the center of the flake. 
As a result the spectrum becomes linear in Eq signifying 
the increased sensitivity to the applied field in a gated 
bilayer. 



II. HYDRODYNAMICS OF A BILAYER 

In this section we describe the formalism and derive 
the plasmon spectrum in a gapless bilayer. In part A the 
equilibrium charge density po(x) is derived within the 
Thomas-Fermi approximation. In part B we obtain an 
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intcgrodiffcrential equation (jSJ for 5p(r, t) and solve it in 
the short-wavelength limit. The solutions have a pseudo- 
continuous spectrum Eqs. (fl~5|) and ([T6|) . which becomes 
discrete, see Eq. (f2U)l and Eq. after logarithmic sin- 
gularities are regularized. The limit of long wavelengths 
is discussed in part C and a surprising minimum in the 
frequency of odd modes is noticed at intermediate wave- 
lengths ~ d. This behavior is further discussed in part 
D. 

We begin with a description of the mean induced 
charge profile po(x) and fluctuation Sp(r,t). The applied 
electric field induces charge density of electrons (or holes) 
Po(x), which has to be found by taking the electric field 
of the induced charges into account self-consistently. The 
induced carrier density is easily estimated by the order of 
magnitude. Consider that the typical distance between 
the charges is denoted by Ie ■ The equilibrium is reached 
when the amount of induced charge becomes sufficient 
to balance the applied external field, e/l E ps Eq, which 
implies for the average density, |po| ~ ~ Eq. In the 
case when Ie <C d there are many induced charges across 
the width of the flake and a continuous description can 
be used. In order to apply the semiclassical approxima- 
tion one more condition has to be satisfied. In particular, 
the plasmon wavelength A has to be much greater than 
the Fermi wavelength of the induced carriers. The latter 
can be estimated from k% ~ |po|/e ~ We conclude 

that it is necessary to ensure that 

V 

Under these conditions macroscopic hydrodynamic 
equations of the charged liquid can be applied. For the 
total electric field we have 

EM) = E -V / dVf^. (2) 
J |r — r'| 

The fluctuating density p(r, t) obeys the charge conser- 
vation law 

p(r,t)+V- J(r,f) =0. (3) 

While the two relations @ and ([3]) are quite generic, 
the remaining equation for the electric current is system 
specific. For single-layer graphenc it has been discussed 
in Refs.[l^andU The band structure of bilayer graphene 
consists of four bands that originate from the coupling 
of two Dirac cones of the two layers.— The two outside 
bands arc separated by a large band gap ~ 0.6 eV and 
will therefore be ignored here. The remaining two bands 
touch each other (in the absence of interlayer bias) and 
are parabolic, ±p 2 /2m, see Fig. 1, with the effective mass 
m = 0.05m , where m is the vacuum electron mass. For 
parabolic bands we find (sec Appendix A for derivation) 

j = \A ( eE - Va*) ; (4) 



(a) E 




FIG. 1: (Color online) (a) Bilayer graphene flake of width 
2d is placed in electric field that separates p and n regions. 
Guided plasmons propagate along the y axis and decay in 
the x direction, (b) Schematic picture of the electron band 
structure: in equilibrium the sum of the electrostatic potential 
and kinetic energy of electrons (at the Fermi level) is constant. 
Vertical arrows indicate the possibility of optical absorption 
with a given frequency: while possible near the center of the 
flake (green non-crossed arrow) , such transitions are forbidden 
(red crossed arrows) when both initial and final states are 
empty or occupied. 

here the chemical potential is related to the induced 
charge density via p = — 7rp/2me, where the fourfold 
(spin and valley) degeneracy is taken into account. The 
form of Eq. (j4]) is readily recognized from the usual ac 
Drudc conductivity a = ie\p\/mu) if one notes that the 
expression in the brackets is simply the gradient of the 
electrochemical potential. The main difference comes 
from the fact that the effective "Drude conductivity" can 
depend on the coordinates (and time) via \p\. 

A. Mean charge distribution 

In equilibrium electric current is absent, J = 0, and 
the mean density of charges po(x) obeys the equation 

f d i i i\ X + X 1 7TOR 

E x + 2 dx'p (x')ln- 7 T + ^ i Po(x) = 0. (5) 

Jo \X — X \ z 

This is the Thomas-Fermi equation where the first two 
terms describe electrostatics of an ideal metal strip. The 
last term takes into account that the screening radius 
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is finite, the latter being of the same order as the Bohr 
radius ag = h 2 /me 2 m 10.6 A. When the width of the 
flake is much greater than the Bohr radius, d 3> as, 
which is the case for any practical situation, the last term 
in Eq. ([5]) is negligible and the solution to the remaining 
integral equation is simply 



Po(x) 



EnX 



Vd 2 



(6) 



This solution is known^S from the method of conformal 
mapping for the two-dimensional Laplace equation, but 
it is most easily verified by a direct substitution into 
Eq. ([5]). Expression is applicable everywhere except 
very close to an edge, where in principle the last term in 
Eq. (|5|) would provide rcgularization of the square-root 
singularity in the mean density. Such a procedure, how- 
ever, strictly speaking, would exceed the accuracy of our 
semiclassical treatment. Indeed, as follows from Eq. (|6|), 
the last term in Eq. §5§ becomes comparable with the 
first two at d — x ~ a 2 B /d, which is a length of the order 
of lattice spacing (or even smaller) , where fully quantum- 
mechanical treatment is warranted. Fortunately for our 
purposes, the singularity in Eq. ^ is integrable and its 
presence will not affect the subsequent calculations. 



B. Plasma oscillations 



plasmon oscillations extend over distances of the or- 
der of their wavelength in the transverse x direction, 
the limits of integration in Eq. (jHJ can be extended 
to infinity while the mean density approximated with 
Po{x) = Egx/d. The plasmon momentum can then be 
conveniently scaled away with the substitution qx = £. 
Furthermore, using the equation for the modified Bcsscl 
function, tKq(\t\) + Kq(\t\) = tKq(\t\), it is convenient 
to rewrite the differential operation in Eq. (JSJ as follows: 

- |£l)#o(K - = C'sgn £(3F - - £'!)■ 

Equation jSJ then becomes (£ = qx): 
2 _ 2eE J d 2 

/OO 
derxKWK -£'!) = (), (9) 
-oo 

that in the Fourier representation acquires the form 



co 2 X (k) 



4neE { 



P 



dk' VT 
2^T "fc 7 



k' 2 dx(k') 



dk' 



(10) 



where P stands for the principal value of integral. The 
integral equation (fTU|) should determine the discrete spec- 
trum of plasmon eigenvalues uj„. As expected, from the 
symmetry of the system, the solutions possess definite 
parity. It can be seen that Eq. (fit)]) has even solutions 



Plasmons are charge density oscillations propagating 
on top of the equilibrium profile , Eq. (jB]), 



P = Po(x) + Sp(r,t). 



(7) 



Equations ©-Q can be linearized with respect to the 
density variation 6p(r,t). The latter will be taken in 
the form of a plane wave propagating along the junc- 
tion, 8p(r,t) = x( x ) cxp(icot — iqy). By substituting for 
electric field E and current J in Eq. we obtain the 
following integrodiffcrcntial equation for the oscillating 
density profile, 



2 / \ 2e 

w x\x) H 

m 



^|po(*)|^-<? 2 M*)| 



dx' X (x')K o {q\x-x'\)=0, 

-d 



(8) 



where Kq is the modified Bessel function of the second 
kind. Note that in deriving Eq. (j8]) we neglect the V/i 
term in Eq. Q. As can be readily seen, its contribution 
is small by the parameter as / A -C 1 • 

Equation (JSJ is reminiscent of the equation for single- 
layer graphene^ where \/\po(x)\ takes the place of 
|po(^)|- This difference makes the solution of the bi- 
layer problem both simpler and trickier. In the most 
interesting case of a wide strip, A <C d, the confine- 
ment of plasmons originates from the gradient of the 
equilibrium charge density. Since, as we find below, 



X^a \k) = cosh cos^a arcsinh kj , (11) 

and odd solutions 

X^a \k) = isinh sin^a arcsinh kj, (12) 

where a is an arbitrary positive number. This follows 
from the following relations 



P / dt 



cosh £ sin(crf;) 7rcos(ar) 



P / dt 



sinht — sinhr tanh(7ra/2) 
coshi cos(ai) 



(13) 



sinh t — sinh r 



7rtanh(7ra/2) sin(ar), (14) 



that can be proven by calculating the corresponding in- 
tegrals with the help of adding an infinite semicircle in 
the complex plane and summing over the residues; see 
Appendix B. The corresponding spectrum of eigenvalues 
is gapped for the even modes, 



"1(a) 



2TreE, 



o 



md tanh(-7ra/2)' 

and gapless for the odd modes, 

2 . . 2ireEo 

cj_(a) = — atanh(7ra/ 2). 

md 



(15) 



(16) 
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Finally, we note that the even and odd solutions obey a 
very simple relation in the real space, 



X i-)(C) = sgn(Oxi +) (e)- 



(17) 



This can be obtained by noticing that X^{k) and 
X (k) are related by the Hilbert transform: the direct 



substitution of the solutions ([TT|) and (fT2j) into Eq. (fit) 
yields 



Xi +) (k) 



-iP 



dk' xt\k') 
7r k' — k 



(18) 



which is the Fourier transform of Eq. (|17[) . 

The most surprising feature of the obtained solutions 
is the continuity of the spectrum of plasmon frequencies, 
which is in a seeming contradiction to the fact that plas- 
mon modes are localized in the transverse direction. To 
elucidate the physical origin of this finding let us find the 
explicit form of our solutions in real space. Due to Eq. 
(fT7|) it is sufficient to consider Xqs (£) f° r positive argu- 
ments, £ > 0. As shown in Appendix B the solution is 
given by the modified Bcsscl function of the imaginary 
order, 



f , c^r^u-t , \ asinh(7ra) 
x / dte'^ CObht cos(at) = ^ — L K ia (£). (19) 



2< 



For large distances, £ ^> 1, the asymptotic behavior is 
exponential, 



Xa (0 K £3/2 ' 



(20) 



We indeed obtain that guided plasmons are localized on 
the distances of order of their wavelengths ~ A. 

The behavior at small distances, £ <?C 1, is more pecu- 
liar: 



my 



2£ r(l + ia) 



(21) 



We observe that at small distances the solutions Eq. (TT9"]) 
have infinitely many nodes, which is formally responsi- 
ble for the continuity of their spectrum. Noteworthy is 

the similarity between our charge density and the 
wave function of a quantum-mechanical particle moving 
in the attractive potential V(£) = — a 2 /£ 2 in two dimen- 
sions, the situation that leads to a particle falling on the 
attraction center^ It is therefore obvious that oscilla- 
tory behavior at £ — > is an artifact of the semiclas- 
sical approximation. The latter fails at small distances 
comparable with the Fermi wavelength (which gets large 
closer to the p-n junction line). Thus such fast oscilla- 
tions are unphysical and would not have appeared in the 



fully quantum-mechanical treatment of electrons in elec- 
tric field. Fortunately, the oscillatory behavior is only 
logarithmic and can be easily regularized. This can be 
performed by noting that the solutions Xa(£) smoothed 
over these fast unphysical oscillations vanish at £ — > 0. 
The vanishing is most easily seen from the density accu- 
mulated at small distances: using Eq. ([2"TT) we find 



dexi +) (0=A cos [a In (£/2) - /?] , (22) 



where /3 is the phase of the complex quantity T(l + ia). 
The accumulated charge therefore oscillates with con- 
stant amplitude A around zero value. We therefore im- 
pose the regularization requirement that the smoothed 
solution should vanish at distances smaller than some a 
(see below) where the semiclassical approach ceases to be 
valid. From Eq. (|T9"| we obtain the equation 



K ia (qa) = 0, 



(23) 



which determines the set of discreet values of a. Note 
that the choice of Eq. (f2"5)) is natural for both the even 
and odd modes as the replacement of the oscillatory func- 
tion Ki a (£ ) with any smoothed nonzero value would have 
led to the nonintcgrablc singularity in Eq. (|T9|) at £ = 0. 

Equation (|23p constitutes the quantization condition 
that together with Eqs. (jlip and (TT2")) determines the 
spectrum of guided plasmon modes for wavelengths 
shorter than the width of the flake. For a « 1 we ob- 
tain in the logarithmic accuracy the following analytic 
approximation: 



In 



2e~i 
qa 



1,2,3, 



(24) 



where 7 = 0.58 is the Euler's constant. The first three 
eigenfunctions X«» (£); n = 1> 2, and 3, are plotted in 

Fig. m 

The cutoff distance a can be estimated as follows. The 
semiclassics fail at the distances of the order of the Fermi 
wavelength, a ~ kp 1 . The latter, however, is a function 
of the coordinate, kp ~ (Egx/de) 1 / 2 , cf. Eq. ©, and 
should be taken at the same distance, x ~ a. We there- 
fore obtain that a ~ ^(d/^E) 1 ^ 3 , i.e., that the cutoff is 
mostly given by the electric length defined in Eq. (TTJ and 
exceeds the latter only by virtue of the factor (d/ls) 1 ^ 3 - 



C. Long wavelengths, q <g; 1/d 

When the plasmon wavelength becomes comparable 
with d (or exceeds it) the oscillating electric field extends 
beyond the width of the flake. The integral equation 
(|8l) has to be solved with the explicit boundary condi- 
tion (unimportant previously) requiring that there is no 
particle flow across the edges of the system, 



J x (±d) = 0. 



(25) 
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FIG. 2: (Color online) The eigenfunctions x&t'CO. Eq. (fT9)l , 
for three lowest values n = 1, 2, and 3 and ga = 0.2. The 
first mode (red) is node-less while the second one (blue) has 
just one node and the third (green) has two nodes and the 
sharpest peak. For the sake of convenience the functions 
2\a (£)/ cosh (-Ka/2) are plotted. The corresponding eigen- 
values for a\ = 1.44, Q2 = 2.49, as — 3.42 are found from 
Eq. (Tig). 



In the limit q — > the properties of plasmon spectrum 
can be elucidated without the exact solution. Eq. © 
becomes 



^ 2 (o)x(C) 



2eE d 



ICI 



md d( y/i - 



, X(C) 

c-c 



(26) 



where £ = x/d. This equation has one zero eigenvalue, 
wi(0) = 0, which corresponds to the function 



x(C) 



(27) 



The solution (|27[) simply describes charge distribution in 
the equipotentially charged metallic strip. For finite but 
small q this mode develops into the usual one-dimensional 
plasmon with the spectrum (cf. Ref. H for a similar- 
discussion of the single-layer case) 



W?.(g) = Cq 2 ln(l/qd). 



(28) 



Here we denote by cu n ± the eigenfrequency of an 
even/odd mode with n—1 nodes across the half width of 
the flake (0,d). In the short- wavelength limit it simply 
corresponds to u>±(a n ); cf. Eqs. (|15p and (|T6"|) . The con- 
stant C is proportional to the total number of induced 
charges per unit length along the y direction and can 
be most simply found by integrating the original Eq. ([5)1 
across the width of the graphene strip. Using the ap- 
proximation Ko(q\x — x'\) = — ]nq\x — x'\ and notic- 
ing that the term containing the total derivative van- 
ishes by virtue of the boundary condition (|23|) , we obtain 
C = AeE^d/m. The logarithm in Eq. (|2"5)) originates from 
the long-range nature of Coulomb interaction. 



Naively, one would expect the lowest-lying odd plas- 
mon solution to be (1—), i.e., to change its sign once, at 
the center of the flake, C, = 0, but to keep the same sign 
across the half width of the sample. Remarkably, Eq. (|2"rJ)) 
does not admit an odd solution without at least one zero 
in the domain (0, d), which also obeys the boundary con- 
dition Eq. (|2"5|) . To see this, let us integrate Eq. (|2l))) over 
£ from to 1. The boundary condition (|2"5j) implies that, 
as C — > 1) we get 



■4(1) oc 



, X(C') 



c-c 



= o. 



(29) 



C-s-i 



Then we arrive at the relation 
w 2 (0)md r ... |C| 



2eE n 



dCx(C) = lim 



, x(C) 



(30) 



In physical terms this means that in order to have a net 
charge across the half width, J* Q d(x(C) 0, there should 
be a nonvanishing current across the p-n junction, given 
by the right-hand side in Eq. (|3"0"|) . Nonvanishing of the 
limit in Eq. (|30[) implies in its turn that at the junction, 
C = 0, the induced electric field must be singularly strong, 



E(C -> 0) cx / d( 



, x(C) 

c-c 



(31) 



To create such a strong field, the plasmon fluctuation 
should have a J- function singularity at £ = 0, x(C) (x 
5(C)- Such a situation is unphysical and in any case in 
conflict with the assumption that x(C) is an °dd function. 
Thus we conclude that j Q d£x(£) = 0, and the lowest- 
lying odd plasmon eigenmode (2—) should have at least 
three nodes across the width of the sample rather than a 
single node as might be intuitively expected. 

The mode (f2"B")) is gapless because its electric poten- 
tial is uniform across the flake. All other modes have 
nodes and therefore their frequencies do not tend to zero 
in the limit q — s- but instead have energy gaps ui n (0) 
determined by Eq. (|2f)|) . Without explicitly solving it we 
conclude from scaling that for all plasmons except the 
lowest even mode, 



LO„(q 0) = const(n) x 



eE 



(32) 



In other words the plasmon spectrum at large wave- 
lengths is determined by the same energy scale eEo/md 
as in the case of short wavelengths (fl~5j) and (fT6|) . 



D. 



Mode order reversal at the intermediate 
wavelengths, q ~ 1/d. 



The results obtained so far have shown surprising re- 
versal in the order of even and odd solutions. At Ions; 
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FIG. 3: (Color online)The sketch of even plasmon (red mono- 
tonic curves) and odd plasmon (blue curves with minima 
around 1/d) frequencies for n = 1 (solid lines) and n = 2 
(dashed lines). The three regions, q <C 1/d, q ~ 1/d, and 
1/d <C q < 1/a are denoted by (a), (b), and (c), respectively. 
With decreasing the wavelength through A ~ d the nth odd 
mode at q <C 1/d develops into the n — 1th odd mode at 
q ~ 1/d, e.g., (2—) becomes (1 — ). 



wavelengths, A 3> d, shown as domain (a) in Fig. [31 
the lowest energy mode is the gapless plasmon wi(g), 
Eq. (|28[), However, in the short- wavelength limit, A <C d, 
domain (c), the order is reversed as w+(a„) > w_(a„), 
cf. Eqs. (fT5j) and ([TB")) . Formally, as the short- wavelength 
cutoff is made smaller, a — ¥ 0, the values of the lowest a n 
decrease and the odd frequencies w_(cv n ) approach zero; 
cf. Eq. ([16)) . Such a change in the order of even and odd 
modes can be explained with simple physical picture of 
charge distribution illustrated in Fig. @] for the first sev- 
eral modes. The lowest even mode's charge distribution 
is nodeless across the width of the flake, and in the case 
of long wavelength A ^> d, Fig. Eta), produces a weak 
restoring longitudinal field E oc 1/A 2 , resulting in the 
gapless spectrum ([2"5)l . At shorter wavelengths the inter- 
mediate domain (b) A ~ d, the distribution of charges 
for the lowest odd mode resembles a checkerboard with 
any plaquette of four charges composing a quadrupolc. 
In the lowest even mode, on the other hand, the plaque- 
tte consists of two uncompensated dipoles that produce 
stronger electric field. We therefore conclude that around 
A ~ d the first odd plasmon becomes the lowest mode of 
the system. 

For intermediate wavelengths A ~ d, the mode (1—) is 
no longer forbidden since J Q d(x(C) 0, as the charge 
transport occurs along the y direction. The (1—) mode 
evolves from the (2—) mode at q = and at A ~ d has the 
frequency below that of the (1+) plasmon. This rever- 
sal happens because the "checkerboard" pattern of the 
(1 — ) mode significantly reduces the electric field com- 
pared with the (1+) arrangement. 

It appears that the reason for the change in the num- 
ber of nodes in odd solutions is the linear decrease of the 
conductivity of a bilayer graphene near the line x — 
that prevents currents from flowing across that line so 



that the zero- node profile (1—) can exist only when the 
wavelength becomes short enough, ~ d, for the longi- 
tudinal currents (along the y axis) to be able to signif- 
icantly affect the charge distribution. In contrast, in a 
monolayer graphene the conductivity vanishes only as 
the square root of x and therefore currents across the 
junction can potentially remain finite if the fluctuating 
field has a 1/ \f\x\ singularity, which is the case for odd 
solutions in a monolayer graphene £ 

Let us emphasize that in plasmonic systems there is 
no inherent reason for even/odd modes to always have 
the same order. In particular, in 3D metal films the even 
mode has lower frequency— while in a system of two 
2D electron layers separated by a dielectric the situation 
is reversed^ What makes guided plasmons in graphene 
bilayer peculiar is the crossover between these two cases 
for different wavelengths. 



E. Comparison with the case of a monolayer, Ref. 8 

It is useful to compare our findings for a bilayer 
graphene with the results of Ref. [1 for a monolayer. In 
the latter case at short wavelengths q ^ 1/d the solu- 
tions were doubly degenerate with a pair of even and odd 
modes having the same frequency. The frequencies were 
oc q 1 / 4 while in the bilayer case they are [see Eqs. ([T5j) and 
([TBI ] virtually independent of q [up to the weak logarith- 
mic dependence in Eq. (|24p j. In addition, the dependence 
on the strength of the electric field creating the p-n junc- 



tion was oc Eq^ 4 compared with the Eq ,z dependence for 
a bilayer. 

The above-mentioned monolayer degeneracy was lifted 
in the long- wavelength limit q < 1/d though not without 
some peculiarities. In particular, the order of the first 
four modes was < ui\- < 0J2- < <^2+- In bilayer 
graphene the situation is more dramatic. From Eq. (|16[) 
we can see that, provided that hil/qa 1, the frequen- 
cies of many odd modes decrease (soften) considerably 
with increasing the wavelength (decreasing q) in the re- 
gion of intermediate wavelengths q ~ 1/d. However, that 
decrease in frequency does not persist at q — > as the 
odd modes remain gapped there. The second unexpected 
finding is that the number of nodes in the odd modes in- 
creases by 1 with increasing the wavelength. This feature 
originates from the fast suppression of the conductivity 
of the bilayer near x = essentially decoupling electric 
currents in the p and n domains of the flake. 
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III. EXPERIMENTAL IMPLICATIONS 

The main obstacle in experimental observation of low- 
dimensional plasmons is that the latter are typically gap- 
less. This makes it impossible to simply convert a pho- 
ton into a plasmon with the conservation of both energy 
and momentum, so that more complicated experimen- 
tal geometries are necessary^ In the case considered in 
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ric modes are likely to be excited. In contrast, symmet- 
ric modes have electric fields that are odd and therefore 
they do not couple to the radiation polarized along the 
x axis. This can be seen in Fig. 2] where the electric 
field of the (2+) mode along the x axis in one half of the 
flake is the mirror reflection (along x = 0) of the elec- 
tric field in other half. In principle, symmetric modes 
could be excited by the longitudinal (along the y axis) 
polarization. However, due to the conservation of mo- 
mentum along the y axis and small values of momentum 
of infrared photons, the phase space for such processes 
is rather restricted. Such absorption has characteristic 
frequencies A = yf AeE$ / md. Using the parameters of 
bilayer graphene we find 
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FIG. 4: (Color online) Schematic charge distribution for the 
lowest plasmonic modes corresponding to different regions in 
Fig- E (a) For long wavelengths, A > d, the lowest mode 
(1+) is the quasi-one-dimensional plasmon with symmetric 
nodeless (across the flake) charge distribution and electric 
field pointed mostly along the y direction. The magnitude of 
electric field vanishes with increasing the wavelength, in agree- 
ment with the gapless spectrum; Eq. (|28p . The second mode 
(2—) is odd and gapped and has three nodes. A mode with a 
single node (1 — ) is forbidden by the condition J^dxx( x ) = 0, 
valid for any gapped mode at q — s> 0. The origin of that con- 
dition is the suppression of charge transport across the line 
separating p and n regions. The third mode (2+) is even and 
has two nodes. The electric field in the modes (2—) and (2+) 
is mostly pointed along the x direction. The field is stronger 
for the mode (2+) (this fact could be easily inferred qualita- 
tively from the picture of charge distribution), thus leading 
to tJ2- < 0J2+- (b) For intermediate wavelengths A ~ d, the 
mode (f — ) is no longer forbidden and the (1 — ) mode at A ~ d 
evolves from the (2—) mode at q — and at A ~ d has the 
frequency below that of the (1+) plasmon. This reversal hap- 
pens because the checkerboard pattern of the (1—) mode has 
lower electric field than the (1+) mode. 



the present paper, plasmons confined near graphene bi- 
layer p-n junctions are in general gapped, cf. Eqs. (|15|) 
and (|32p. and therefore straightforward absorption of 
long-wavelength (infrared) radiation is possible. Since 
the wavelength of the infrared radiation of interest to 
us reaches millimeter range the corresponding electric 
field is virtually homogeneous, so that only antisymmet- 



A w 2.7 eV x 




(33) 



For electric fields, ~ 10 6 V/m, the value of electric length 
is Ie ~ 40nm. Assuming d ps for the size of the 

sample we obtain from Eq. ([33"]) that A ~ 2 meV. We 
therefore expect a threshold signature in the absorption 
spectrum of infrared radiation at these frequencies, which 
are very sensitive to the magnitude of the electric field 
creating the p-n junction. 

We emphasize that the plasmon absorption happens 
in addition to the intcrsubband electron-hole absorption. 
The latter, however, is a smooth function of frequency. 
Indeed, interband absorption of a photon with frequency 
oj is possible only as long as the chemical potential (the 
distance to the Fermi level from the band degeneracy 
point p = 0) does not exceed hu/2; see Fig. 1(b). 
The former, expressed via the charge density, Eq. ((BJ, 
is simply 7r|p|/2me = ir Eq\x\ / 2me\f d 2 — x 2 . Thus the 
corresponding transitions are allowed within a strip of 
| a; | < dj s/T+ (n Eq / meui) 2 . The intensity of the absorp- 
tion is independent of frequency and simply determined 
by the intersubband ac ( "minimal" ) conductivity^ of bi- 
layer, e 2 /2h. The frequency dependence of single-particle 
absorption is therefore simply determined by the width 
of the absorbing region, 



A e - h (uj) oc 



(34) 



This expression describes a smooth background that ex- 
ists in addition to the usual intrasubband Drudc absorp- 
tion. Its characteristic frequency Ti 2 /ml 2 E depends lin- 
early on the applied electric field and should be easily 
distinguishable by varying the electric field from the plas- 
mon contribution whose frequency Q33D is proportional to 
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E, 



IV. GATED BILAYER 

The advantages of bilayer graphene for applications 
lie in the possibility of inducing a bandgap with 



the interlayer bias, studied both theoretically^— and 
experimentally^ - — Since this is typically performed via 
the nearby metallic gates, we now discuss how the prop- 
erties of guided plasmons are modified by the presence 
of such a gate, which is assumed to be positioned a dis- 
tance D from the bilayer. In addition, the presence of 
a bandgap leads to the appearance of a neutral strip of 
width 2h; see Fig. [SJ In this section after explaining 
how the general equations are changed, we consider three 
cases: (i) short wavelengths, q^ 1 «D, (ii) intermediate 
wavelengths, D <C q^ 1 <C h, and (iii) long wavelengths, 
/i< q^ 1 < d. 

Two modifications have to be made to our hydrody- 
namic equations. In particular, the induced electric field 
includes the additional contribution from image charges, 



E(r,*) = E -V / d 2 r'p(r',t) 



AD 2 



(35) 



The second modification should incorporate the pres- 
ence of a gap 2U in the energy spectrum, e(p) = 
±\/[/ 2 + {p 2 /2m) 2 . The relation between the local 
chemical potential and the induced charge density now 
becomes 



-sgn (p) \ U 2 + 



\2me 



) 



(36) 



The existence of a gap in the spectrum means that 
separation of electrons and holes is possible only when 
the applied in-plane electric field exceeds some value, 
Eq > U /erf. Even above this threshold p and n regions 
are separated by a neutral strip of width 2h = 2U/eE, 
as shown in Fig. 2. The mean distribution of the in- 
duced charges po(x) is determined by the Thomas-Fermi 
equation 



E x + sgn(x) 



B „2 



dx' po(x') In ■ 



(37) 



The second term under the square root can (similarly to 
Sec. II) be neglected everywhere except very close to the 
flake's edges [more precisely, as long as cibPo(x) is less 
than either Eqx or U/e}. The resulting linear equation 
has been solved in Ref. [2Jj in relation to the polarization 
of a nanotube array in the external electric field, yielding 



Po(x) = sgn(x) E 



^B(\x\-h), (38) 



x 



The equation of motion for plasmon oscillations is also 
modified by the gap in the spectrum, as explained in 
Appendix A; see Eq. (|A5|) . Solution of the corresponding 
general equation for the oscillating density is beyond the 




FIG. 5: (Color online) Electron band structure across the bi- 
layer with the energy gap of 2U. The p (blue line starting 
from the point where the second lower peak on the left side 
touches the Fermi level and ending at the left end of the flake) 
and n regions (red line on the right side of the axis) are sep- 
arated by a neutral (TV) strip (black double-headed arrow) of 
width 2h = 211/ eE determined by the energy gap and the 
slope of the scalar potential. 



scope of the present paper. However, great simplification 
occurs if the gap is not very small, so that the width of 
the neutral strip is much greater than the Bohr radius, 
h 3> as- In this case the second term in the denominator 
of Eq. (|A5|) is negligible compared with U, which amounts 
to the substitution 



I Pol 



2meU 



(39) 



in Eq. ([5]). The second change to Eq. comes from tak- 
ing into account the image charges from Eq. (|35[) . This 
is done by replacing Ka(q\x — x'\) with 



K (q\x -x'\)- K [q^{x-x') 2 +4D 2 } 



dk 



Below we elucidate plasmon dispersion for the case of a 
wide flake and a gate situated in close proximity to it, 
L> « h < d. 

(i) Short wavelengths, q^ 1 In this case the effect 

of the image charges is negligible. Plasmons propagate 
along a boundary between a charged (p or n) domain and 
a neutral strip (TV) and decay over a distance short com- 
pared to h. Thus modes localized near the p-N boundary 
are exponentially weakly coupled to n-N modes. The cor- 
responding plasmon frequencies, 



E 2 h 
m 2 d 2 U ' 



(41) 



are therefore doubly degenerate and independent of the 
wavelength. 

(ii) Intermediate wavelengths, D <C q^ 1 <C h. The 
wavelength is still short enough to ensure that plasmons 
propagating near the two edges of the neutral strip do not 
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"talk," but the Coulomb interaction is screened by the 
gate. The expression in Eq. f|40[) can be approximated 
with 2itD5(x — x'). Since only small x — x — h <C h 
are of interest, the density in Eqs. (j3"7| and can be 
approximated with (x > 0) 



Po( x ) 



IhEl 



(42) 



By virtue of the 8 function the equation for the oscillating 
charge density becomes a differential one, 



1 d 



x dx V dx 



( 



\C 2 x 



X (x) = 0, 



(43) 



where C 2 



4tt 2 ElhD 



This equation is identical to that 

2 



d 2 m 2 U 

of a two-dimensional "hydrogen atom" with — q 2 playing 
the role of the energy and ui 2 /C 2 the role of the interac- 
tion constant. The corresponding quantization condition 
is simply q — ui 2 /2C 2 (n+l/2), which yields the spectrum 



u) 2 n = (n + 1/2) 



8ir 2 E 2 hD 
d 2 m 2 U 



n = 0,1,2, 



(44) 



The corresponding eigcnfunctions are given by the con- 
fluent hypergeometric function of the first kind, 



Xn(x) 



iFi(-n,l,2 9 5). 



(45) 



Remarkably, the solutions (|45|) do not vanish at x = 
x — h = 0. This indicates that the shape of the boundary 
of the neutral strip, created by "soft" electrostatic con- 
finement, does not remain linear and itself fluctuates due 
to plasmon excitations. 

(iii) For yet longer wavelengths, h <C q^ 1 <C d, the 
presence of the neutral strip is irrelevant as the field of 
plasmons extends over much larger distances. Assuming 
that the wavelength is still much smaller than the width 
d, we can use p 2 = E^x 2 jd 2 '', cf. Eq. ©. The 5 function 
approximation from the previous paragraph still holds, 
yielding 



2<h_ 
di 



2f2 



-i )x = o, 



(46) 



where £ 



qx, and C\ = ^af ■ With the help of the 



substitution, x(Q = f(0/\Z\> Ec i- © is reduced to 



-^+l7^272 -l)/ = 



Cti 



(47) 



which is the one-dimensional Schrodinger equation for 
the potential cx — l/^ 2 . The singularity at £ = needs 
to be regularized by cutting off at £ « qa. Using the 
Bohr-Sommcrfeld condition, 



c 2 e 



1 = 2tt \n 



(48) 



we find the spectrum to be 



tt 4 E 2 D 



2m 2 d 2 U\n 2 (qa) 



(n + 1/2) 2 , n = 0,l,2, 



(49) 



The result (|49p becomes more accurate for n 1, since 
the WKB applicability condition for Eq. (|47| requires 
that uj C\. Still, even for n ~ 1 the expression (j^)) 
captures the correct dependence on the gap U and the 
strength of the external field Eq. 

To conclude the discussion of the gated bilayer, let 
us note that only in the intermediate case (ii) does the 
plasmon frequency demonstrate dependence on the wave 
number, namely when the wavelength is long enough to 
ensure screening by image charges, but still sufficiently 
short compared with the width of the neutral strip. In 
all three regimes the dependence of plasmon spectra on 
electric field E is linear in contrast to a gapless case, 
where the corresponding dependence is square root; see 
Eqs. (HD, CEU), and ([32]). 



qa 



V. SUMMARY 

Bilayer graphene is a gapless (or weakly gapped) sys- 
tem that is charge neutral when undoped. Separation of 
charges, however, occurs with even weak external elec- 
tric fields applied along the plane of bilayer; see Fig. 1. 
The induced charge density Eq. (JB]) follows from the solu- 
tion of the electrostatic problem. The gradient of charge 
density near the p-n junction line x = leads to the 
confinement of charge oscillations (plasmons) resulting 
in their one-dimensional propagation along the junction. 
Physically, confinement in the transverse direction can 
be understood as follows. Since plasmon "stiffness" in- 
creases with increasing the local charge density, plasmons 
of lower frequency tend to be located closer to the p-n 
junction and to decay into the region of higher density. 

Within a continuous hydrodynamic approach for 
charge-density oscillations the plasmon modes are de- 
termined from the integrodiffercntial eigenvalue problem, 
Eq. (jTU)) . The latter allows for exact solutions, Eqs. (fTTj) 
and (|T2"j) , with the discretization of the spectrum being a 
consequence of the rcgularization of the logarithmic sin- 
gularity, Eq. (|2"3")h For small wave vectors q the lowest 
mode is an even solution that is a conventional ID plas- 
mon with logarithmic velocity. At q ~ l/d the order 
is reversed and the first odd solution becomes the mode 
with the lowest energy. 

From the practical standpoint guided plasmon modes 
could potentially be used for "plasmon transistors."— The 
immediate experimental signatures of the guided plas- 
mons, however, can be most directly measured in the 
optical absorption. While the electron-hole intersub- 
band excitations lead to a smooth continuum, Eq. (|34[) . 
the gapped plasmon spectrum of antisymmetric modes 
should result in a threshold behavior at a frequency 
~ eEo I md sensitive to the applied electric field. 
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vp = d/i/dpF and the density of states, v{p) — 2pf/ttvf, 
we obtain 

_ „2 

: (eE-V/i). (A5) 



2m 2 e 



\2me/ 



In the case where there is no gap, {7 = 0, Eq. (HJ) is 
recovered. 



Appendix A: Hydrodynamic equations 

Equations (|3|) and (|4|) can be derived from the Boltz- 
mann equation for the electron distribution function / p , 



at 



v • V/ p + eE ■ 



dfe 
dp 



lift 



pjj 



(Al) 



where the right-hand side is the collision integral^ Since 
the latter conserves the number of electrons, integrating 
the equation over the entire momentum space yields the 
continuity equation ©. Hydrodynamic approximation 
is applicable when a distribution function deviates only 
slightly from a local equilibrium distribution with the 
chemical potential p(r, t) (which is equivalent to specify- 
ing the local density p) . The deviation is due to the drift 
of particles with the average velocity being considerably 
smaller than the Fermi velocity vf , and characterized by 
a local current density J(r,t): 



/ p (r,t) = e(/i(r,t)-e p )+2- 



%(r,t)-e p ], (A2) 



where v{p) is the density of states at the Fermi level. 
Such an ansatz is nothing but the expansion over angular 
harmonics truncated after the first term. Multiplying 
Eq. (|A1|) by ev and performing the same operation, we 
can obtain the equation for electric current, 



§ + e£v(v.V/ p ) 



dp 



(A3) 



The last term represents collision relaxation of the elec- 
tric current and for small currents (linear response) may 
be written as — J/r, where r is the transport mean free 
time. Substituting Eq. (|A2|) into Eq. ([S5]l we get 



- = o ew F^(M)(eE- 

T I 



at 



(A4) 



Appendix B: Integral relations 



Integral relations, Eqs. (jl3[) and (jT4|) . can be estab- 
lished by considering the integral 



P / dt 



sinh t — sinh t 



(Bl) 



This integral can be evaluated by noticing that in the 
upper plane of the complex variable t, the integrand is 
exponentially suppressed on the infinite semicircle, t = 
Re iv , < (p < it, R -> oo; sec Fig. 6(a). Then the 
integral equals the sum of residues at t = r + 27rni and 
t = — t + 7r(2n — where n = 1,2, . . ., times 27ri, plus 
the residue ait — t times ni, which is due to the principal 
integration. After the integral of Eq. (|B1[) is evaluated, 
Eqs. (fT3"|) and (|14p hold as being the real and imaginary 
parts of it, respectively. 

Upon replacing k i— > sinhr and k' i— > sinhi, the inte- 
grodiffcrcntial equation (|10[) acquires the form 



w 2 x(r) = 



AneE 
md 



coshi 



dm 



2ir sinh t — sinh r 



(B2) 



where = x(sinht). Equations (fT3|) and ([M]) imply 
that x(0 = cos(at) and = sin(crf) solve the inte- 
grodiffcrential equation, Eq. (|B2[) . yielding the spectra 
Eqs. (|T5|) and (fT6|) . respectively. Going back to the orig- 
inal variable, t = arcsinhfc, we recover the even and odd 
solutions, Eqs. (fTTj) and (|T2l . which in addition are mul- 
tiplied by factors cosh(7ra/2) and i sinh(7ra/2), respec- 
tively, for the sake of further convenience. This can be 
done without destroying the solutions as these factors are 
constants with respect to the variable k. 

In the remaining portion of this Appendix we establish 
the relation Eq. (fl~9| . Consider the Fourier transform 



2tt j 



X ( a +) (k). 



(B3) 



This expression generalizes the usual Drude conductiv- 
ity onto the case of coordinate- and time-dependent den- 
sity of 2D electron gas. For frequencies exceeding the 
collision rate the second term in the left-hand side can 
be neglected. If the gap is present the spectrum is 
e(p) = \/U 2 + (p 2 /2m) 2 . Calculating the Fermi velocity 



After performing a change, k = sinht, and integrating 
by parts, we rewrite Eq. (|B3j) in the form 



_ acosh(7ra/2) 
Ami; 



oo 
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The integrand in Eq. (|B4[) is a holomorphic function of 
the complex variable t, so that the integration contour 
can be deformed from real axis to the horizontal line, 
— oo + in/2 < t < oo + «7r/2; see Fig. 6(b). In doing this 
we also notice that the integrals over the vertical parts, 
— oo < t < — oo + iir/2 and oo + in/2 < t < oo, are 
suppressed. Then Eq. (|B4|) turns into the relation 

oo 

2mi J 

— OO 

+«cos(o!i)sinh(7ra/2)]. (B5) 

The first term in the rectangular brackets turns to zero 
as being an odd function of t, and we finally arrive at the 
integral representation, 



xi +) (6 



a sinh(7ra) 
2^1 



die 



-£ cosh t 



cos(at). (B6) 




b) 



in/2 



FIG. 6: (a) Contour in the complex plane for the calculation 
of the integral (|B1[) ; the singularities in the integrand are 
indicated by the dots; (b) deformation of the contour from 
the real axis for the calculation of the integral (|B4[) . 
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